# Load packages
```{r}
library(rethinking)
```

# Load observational data and simulation data from model 3
```{r}
load("./data/d_sol_cw_w_in.RData")
load("./data/d_sim_in.RData")
load("./data/model_out_mc.RDATA")
load("./data/modeldata_out_mc.RData")
load("./data/model_out_m3.RDATA")
load("./data/modeldata_out_m3.RData")
```

# Create simulation data for Algae collider
```{r}
## Set parameters
a_A <- 3
b_A_S <- 0.7
b_A_C <- -0.2

d_sim_mc <- d_sim_m3

## Generate simulated collider Algae variable
d_sim_mc$A <- rnorm(length(d_sim_mc$Time), a_A + b_A_S*d_sim_mc$Sol_std + b_A_C*d_sim_mc$C_MW, 0.4)
d_sim_mc$A_std <- standardize(d_sim_mc$A)

windows(); plot(d_sim_mc$Sol)
windows(); plot(d_sim_mc$C_MW)
windows(); plot(d_sim_mc$A)

save(d_sim_mc, file = "./data/d_sim_mc_in.RData")
```

## 3. Define model m3_A: T_S -> C_MW <- Sol with interaction T_S -> Sol
```{r}
# Define the input data
load("./data/d_sim_mc_in.RData")
d_in <- d_sim_mc

dlist <- list(
    C_MW = (d_in$C_MW),
    A_std = (d_in$A_std),
    Sol_std = (d_in$Sol_std),
    T_S_std = (d_in$T_S_std)
)

model_mc <- ulam(log_lik=TRUE,
    alist(
    # Sol -> C_MW <- T_S model and association C_MW - Algae
    C_MW ~ dnorm(mu, sigma),
    mu <- a_C + b_C_S * Sol_std + b_C_T * T_S_std + b_A_C * A_std,
    a_C ~ dnorm(0, 1),
    b_C_S ~ dnorm(0.5, 0.5),
    b_C_T ~ dnorm(0.5, 0.5),
    b_A_C ~ dnorm(0.5, 0.5),
    sigma ~ dexp(1),

    # Sol -> T_S model
    T_S_std ~ dnorm(nu, tau),
    nu <- a_T + b_T_S * Sol_std,
    a_T ~ dnorm(0.5, 1),
    b_T_S ~ dnorm(0.5, 0.5),
    tau ~ dexp(1),

    # Sol -> Algae model
    A_std ~ dnorm(xi, phi),
    xi <- a_A + b_A_S * Sol_std,
    a_A ~ dnorm(2, 1),
    b_A_S ~ dnorm(0.5, 0.5),
    phi ~ dexp(1)
    ), data = dlist, chains = 2, cores = 2
)

precis(model_mc, prob=0.9, digits = 3)
```

# Save the data
```{r}
save(model_mc, file = "./data/model_mc_out.RDATA")

data_mc <- extract.samples(model_mc)

save(data_mc, file = "./data/modeldata_mc_out.RData")
```

